home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / RNG.ZIP / RAN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-19  |  14.8 KB  |  474 lines
  1. /* ran.c - functions to initialize random number generator and random number
  2. ** generator itself.
  3. **
  4. ** Module entry points
  5. ** -------------------
  6. **
  7. ** setup()  ----> N.B.  MUST be called before irand() or poiss()
  8. ** ran()    ----> macro defined in ran.h
  9. ** irand()
  10. ** poiss()
  11. ** binom()
  12. **
  13. **
  14. ** Either RAN_KNU or RAN_MAR must be defined before compilation, but
  15. ** only one can be defined.  If RAN_KNU is defined, a subtractive
  16. ** generator from volume two of Knuth is used.  If RAN_MAR is defined,
  17. ** a more complex generator discovered by Marsaglia is used.  Knuth's
  18. ** generator is about 50% faster than Marsaglia's, but its sequence is
  19. ** not as random.  For many purposes, Knuth's generator is good enough.
  20. ** It is much better than the generators normally supplied with C compilers,
  21. ** and it is often faster.
  22. **
  23. ** setup() must be called before irand() is called.  irand() does
  24. ** NOT check to see that the array it uses has been initialized.
  25. **
  26. ** irand() returns a uniformly distributed integer in [0..NO_NUMBERS)
  27. ** when NO_ZERO is not defined and in (0..NO_NUMBERS) when NO_ZERO is
  28. ** defined.
  29. **
  30. ** The macro ran() (defined in ran.h) returns a uniformly distributed
  31. ** random number in the range [0,1) when NO_ZERO is not defined and in
  32. ** (0,1) when NO_ZERO is defined.  The resolution is 1/CONV_REAL.
  33. **
  34. ** Strictly speaking poiss() should only be called when NO_ZERO is not
  35. ** defined.  But the error will be very small, since it will have an
  36. ** effect on the results in only 1 out of NO_NUMBERS calls, on average.
  37. **
  38. **
  39. ** References:
  40. **
  41. ** Knuth, D. E.  1981.  The Art of Computer Programming.  Vol. 2,
  42. **    Seminumerical Algorithms.  Addison-Wesley, Reading Mass.
  43. **
  44. **
  45. ** Author:  Kent E. Holsinger
  46. **
  47. ** Revision history
  48. ** ----------------
  49. **
  50. ** Version 1.0 -- 24 January 1990
  51. ** Original version of Knuth additive random number generator
  52. **
  53. ** Version 1.1 -- 24 January 1990
  54. ** ran() changed to macro operating on long
  55. **
  56. ** Version 2.1 -- 31 December 1990
  57. ** Added Marsaglia generator with conditional compilation
  58. **
  59. ** Version 2.3 -- 18 December 1991
  60. ** Converted both generators to use unsigned longs
  61. ** Editorial changes to internal comments
  62. **
  63. **
  64. */
  65.  
  66. #include <stdio.h>
  67. #include <math.h>
  68. #include "ran.h"
  69.  
  70.  
  71.                            /* constants used in random number generator */
  72. #ifdef RAN_KNU
  73. #define D      21
  74. #define LASTDIM      55
  75. #define FIRSTDIM  24
  76. #endif
  77. #ifdef RAN_MAR
  78. #define CD     7654321
  79. #define CM     16777213
  80. #define LASTDIM   98
  81. #define FIRSTDIM  34
  82. #endif
  83.  
  84.                            /*
  85.                            ** array of integers used in random number
  86.                            ** generator
  87.                            */
  88. typedef unsigned long INTARRAY[LASTDIM];
  89. typedef unsigned long *POINTER;
  90.  
  91.                            /* file used to store random number seed */
  92. FILE *randfile;
  93.  
  94.                            /* LOCAL VARIABLE DEFINITIONS */
  95. static INTARRAY x;         /* array used to generate random nos. */
  96. static POINTER j_ptr, k_ptr;  /* pointers into array */
  97.  
  98. #ifdef RAN_MAR
  99. static long c;
  100. #endif
  101. static char MOD_ID[] = "%W% %D%";
  102.  
  103.  
  104.                            /* LOCAL FUNCTION PROTOTYPES */
  105. #ifdef __STDC__
  106. static void rknuin(unsigned long m);
  107. static void rmarin(int ij, int kl);
  108. #else
  109. static void rknuin();
  110. static void rmarin();
  111. #endif
  112.  
  113.  
  114.  
  115. /***************************************************************************/
  116. /*                                                                         */
  117. /*                           MODULE ENTRY POINT                            */
  118. /*                                 setup()                                 */
  119. /*                                                                         */
  120. /***************************************************************************/
  121.                            /* Initialize the random number generator  */
  122. #ifdef __STDC__
  123. void setup(unsigned long *seed)
  124. #else
  125. void setup(seed)
  126.    unsigned long *seed;
  127. #endif
  128. {
  129. #ifdef RAN_MAR
  130.    unsigned long temp;     /* to write seed to file for next call */
  131.    int ij, kl;              /* seeds for Marsaglia generator */
  132. #endif
  133.  
  134.    randfile = fopen("RAND.INT", "r");
  135. #ifdef RAN_KNU
  136.    if (!randfile) {
  137.       printf("Enter a random number seed: ");
  138.       scanf("%lu", seed);
  139.    } else {
  140.       fscanf(randfile, "%lu", seed);
  141.       fclose(randfile);
  142.    }
  143.    do {
  144.       if (*seed == 0) {
  145.     fprintf(stderr, "Seed must be greater than zero\n");
  146.     fprintf(stderr, "Enter a random number seed: ");
  147.     scanf("%lu", seed);
  148.       }
  149.    } while (*seed == 0);
  150.    rknuin(*seed);
  151.    randfile = fopen("RAND.INT", "w");
  152.    fprintf(randfile, "%lu", x[LASTDIM - 1]);
  153.    fclose(randfile);
  154. #endif
  155. #ifdef RAN_MAR
  156.    if (!randfile) {
  157.       printf("Enter first random number seed:  ");
  158.       scanf("%d", &ij);
  159.       printf("Enter second random number seed: ");
  160.       scanf("%d", &kl);
  161.    } else {
  162.       fscanf(randfile, "%lu", seed);
  163.       fclose(randfile);
  164.       ij = (int) (*seed) % 31328;
  165.       kl = (int) ((*seed >> 16) % 30081);
  166.    }
  167.    if (ij < 0 || ij > 31328 || kl < 0 || kl > 30081) {
  168.       do {
  169.     fputs("The first random number seed must have a value between 0\
  170.  and 31328.\n", stderr);
  171.     fputs("The second seed must have a value between 0 and 30081.\n",
  172.           stderr);
  173.     printf("Enter first random number seed:  ");
  174.     scanf("%d", &ij);
  175.     printf("Enter second random number seed: ");
  176.     scanf("%d", &kl);
  177.       } while (ij < 0 || ij > 31328 || kl < 0 || kl > 30081);
  178.    }
  179.    rmarin(ij, kl);
  180.    temp = ((x[LASTDIM - 1] >> 16) % 30081) << 16;
  181.    temp += x[LASTDIM - 1] % 31328;
  182.    randfile = fopen("RAND.INT", "w");
  183.    fprintf(randfile, "%lu", temp);
  184.    fclose(randfile);
  185. #endif
  186. }
  187.  
  188.  
  189.  
  190. /***************************************************************************/
  191. /*                                                                         */
  192. /*                           MODULE ENTRY POINT                            */
  193. /*                                 irand()                                 */
  194. /*                                                                         */
  195. /***************************************************************************/
  196.                            /* Returns a random number in [0,NO_NUMBERS) */
  197.                            /* if NO_ZERO left undefined */
  198.                            /* Returns a random number in (0,NO_NUMBERS) */
  199.                            /* when NO_ZERO is defined */
  200. #ifdef RAN_KNU
  201.                            /* Knuth generator */
  202.  /*
  203.   * This is a direct translation of the subtractive generator (p. 171),
  204.   * which is based on the additive generator described on p. 27.  It has
  205.   * a period greater than 2^55 - 1.
  206.   *
  207.   * This implementation takes advantage of the fact that all elements of
  208.   * the array are themselves random numbers to return the element in the
  209.   * array currently pointed to.  Strict implementation of the algorithm
  210.   * would require definition of a temporary variable to hold the return
  211.   * value, which properly is *j_ptr before it is decremented
  212.   *
  213.   * On a Zenith Z-386 each call requires about 13 microseconds with
  214.   * NO_ZERO defined.  Without NO_ZERO defined each call requires about
  215.   * 11 microseconds.  (Results obtained using Microsoft C v6.00a,
  216.   * optimizing with /Ozx
  217.   */
  218. #ifdef __STDC__
  219. unsigned long irand(void)
  220. #else
  221. unsigned long irand()
  222. #endif
  223. {
  224.    *j_ptr -= *k_ptr;
  225.    j_ptr--;
  226.    k_ptr--;
  227.    if (j_ptr < x)
  228.       j_ptr = &x[LASTDIM - 1];
  229.    else if (k_ptr < x)
  230.       k_ptr = &x[LASTDIM - 1];
  231. #ifndef NO_ZERO
  232.    return *j_ptr;    /* simply return *j_ptr for U in [0,1) */
  233. #else
  234.    if (*j_ptr == 0)
  235.       return irand();
  236.    else
  237.       return *j_ptr;
  238. #endif
  239. }
  240. #endif
  241.  
  242. #ifdef RAN_MAR
  243.  /* Marsaglia generator */
  244.  /*
  245.   * C This random number generator originally appeared in "Toward a Universal
  246.   * C Random Number Generator" by George Marsaglia and Arif Zaman. C Florida
  247.   * State University Report: FSU-SCRI-87-50 (1987)
  248.   * C
  249.   * C It was later modified by F. James and published in "A Review of Pseudo-
  250.   * C random Number Generators"
  251.   * C
  252.   * C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
  253.   * C       (However, a newly discovered technique can yield
  254.   * C         a period of 10^600. But that is still in the development stage.)
  255.   * C
  256.   * C It passes ALL of the tests for random number generators and has a period
  257.   * C of 2^144, is completely portable (gives bit identical results on all
  258.   * C machines with at least 24-bit mantissas in the floating point
  259.   * C representation).
  260.   * C
  261.   * C The algorithm is a combination of a Fibonacci sequence (with lags of 97
  262.   * C   and 33, and operation "subtraction plus one, modulo one") and an
  263.   * C   "arithmetic sequence" (using subtraction).
  264.   * C========================================================================
  265.   * This version is based on a C version by Jim Butler, which is itself based
  266.   * on a FORTRAN program posted by David LaSalle of Florida State University.
  267.   *
  268.   * This version takes advantage of the 32-bit long integers in ANSI C to gain a
  269.   * significant speed increase.  On a Zenith Z-386 running DOS 3.2 the
  270.   * floating point version requires about 74 microseconds per call.  This
  271.   * version requires about 19 microseconds per call without NO_ZERO defined.
  272.   * With NO_ZERO defined each call requires about 21 microseconds.  (Results
  273.   * obtained using Microsoft C v6.00a, optimizing with /Ozx.)
  274.   */
  275. #ifdef __STDC__
  276. unsigned long irand(void)
  277. #else
  278. unsigned long irand()
  279. #endif
  280. {
  281.    long uni;
  282.  
  283.    uni = *j_ptr - *k_ptr;
  284.    if (uni < 0) {
  285.       uni += NO_NUMBERS;
  286.    }
  287.    *j_ptr = uni;
  288.    j_ptr--;
  289.    if (j_ptr == x) {
  290.       j_ptr = &x[LASTDIM - 1];
  291.    }
  292.    k_ptr--;
  293.    if (k_ptr == x) {
  294.       k_ptr = &x[LASTDIM - 1];
  295.    }
  296.    c -= CD;
  297.    if (c < 0) {
  298.       c += CM;
  299.    }
  300.    uni -= c;
  301.    if (uni < 0) {
  302.       uni += NO_NUMBERS;
  303.    }
  304. #ifndef NO_ZERO
  305.    return (unsigned long) uni;   /* simply return uni for U in [0,1) */
  306. #else
  307.    if (uni == 0)
  308.       return irand();
  309.    else
  310.       return (unsigned long) uni;
  311. #endif
  312. }
  313. #endif
  314.  
  315.  
  316. /***************************************************************************/
  317. /*                                                                         */
  318. /*                           MODULE ENTRY POINT                            */
  319. /*                                 poiss()                                 */
  320. /*                                                                         */
  321. /***************************************************************************/
  322. #ifdef __STDC__
  323. int poiss(double mean)
  324. #else
  325. int poiss(mean)
  326.    double mean;
  327. #endif
  328. {
  329.    double e_u;
  330.    double u;
  331.    int k;
  332.  
  333.    e_u = exp(-mean);
  334.    u = ran();
  335.    for (k = 1; u > e_u; k++) {
  336.       u *= ran();
  337.    }
  338.    return (k - 1);
  339. }
  340.  
  341.  
  342.  
  343. /***************************************************************************/
  344. /*                                                                         */
  345. /*                           MODULE ENTRY POINT                            */
  346. /*                                 binom()                                 */
  347. /*                                                                         */
  348. /***************************************************************************/
  349. #ifdef __STDC__
  350. int binom(int n, double p)
  351. #else
  352. int binom(n, p)
  353.    int n;
  354.    double p;
  355. #endif
  356. {
  357.    int i;            /* index variable */
  358.    int ct;           /* cumulative number of hits in binomial */
  359.  
  360.    ct = 0;
  361.    for (i = 0; i < n; i++) {
  362.       if (ran() < p) {
  363.     ct++;
  364.       }
  365.    }
  366.    return ct;
  367. }
  368.  
  369.  
  370.  
  371.  
  372.  
  373. /***************************************************************************/
  374. /*                                                                         */
  375. /*                             LOCAL FUNCTIONS                          */
  376. /*                                 rknuin()                                */
  377. /*                                 rmarin()                                */
  378. /*                                                                         */
  379. /***************************************************************************/
  380. #ifdef RAN_KNU
  381.  /* initializer for Knuth generator */
  382.  /* It is based on the one suggested on */
  383.  /* p. 172 */
  384. #ifdef __STDC__
  385. static void rknuin(unsigned long m)
  386. #else
  387. static void rknuin(m)
  388.    unsigned long m;
  389. #endif
  390. {
  391.    double dummy;        /* dummy used to "warm up" generator */
  392.    unsigned long int n;       /* integers used to load generator */
  393.    int i;            /* index variable */
  394.    int item;            /* item storage */
  395.  
  396.    x[LASTDIM - 1] = m;
  397.    n = 1;
  398.    for (i = 1; i < LASTDIM; ++i) {
  399.       item = ((D * i) % LASTDIM) - 1;
  400.       x[item] = n;
  401.       n = m - n;
  402.       if (n < 0)
  403.     n += NO_NUMBERS;
  404.       m = x[item];
  405.    }
  406.    j_ptr = &x[FIRSTDIM - 1];
  407.    k_ptr = &x[LASTDIM - 1];
  408.    for (i = 0; i < LASTDIM; ++i)
  409.       dummy = ran();
  410. }
  411. #endif
  412.  
  413. #ifdef RAN_MAR
  414.  /* initializer for Marsaglia generator */
  415.  /*
  416.   * C This is the initialization routine for the random number generator RANMAR()
  417.   * C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
  418.   * C                                                      0 <= KL <= 30081
  419.   * C The random number sequences created by these two seeds are of sufficient
  420.   * C length to complete an entire calculation with. For example, if several
  421.   * C different groups are working on different parts of the same calculation,
  422.   * C each group could be assigned its own IJ seed. This would leave each group
  423.   * C with 30000 choices for the second seed. That is to say, this random
  424.   * C number generator can create 900 million different subsequences -- with
  425.   * C each subsequence having a length of approximately 10^30.
  426.   * C
  427.   * C Use IJ = 1802 & KL = 9373 to test the random number generator.  The
  428.   * C subroutine RANMAR should be used to generate 20000 random numbers.
  429.   * C Then display the next six random numbers generated multiplied by 4096*4096
  430.   * C If the random number generator is working properly, the random numbers
  431.   * C should be:
  432.   * C           6533892.0  14220222.0  7275067.0
  433.   * C           6172232.0  8354498.0   10633180.0
  434.   */
  435. #ifdef __STDC__
  436. static void rmarin(int ij, int kl)
  437. #else
  438. static void rmarin(ij, kl)
  439.    int ij;
  440.    int kl;
  441. #endif
  442. {
  443.    int i, j, k, l, ii, jj, m;
  444.    float s, t;
  445.  
  446.    i = (ij / 177) % 177 + 2;
  447.    j = ij % 177 + 2;
  448.    k = (kl / 169) % 178 + 1;
  449.    l = kl % 169;
  450.  
  451.    for (ii = 1; ii <= LASTDIM - 1; ii++) {
  452.       s = 0.0;
  453.       t = 0.5;
  454.       for (jj = 1; jj <= 24; jj++) {
  455.     m = (((i * j) % 179) * k) % 179;
  456.     i = j;
  457.     j = k;
  458.     k = m;
  459.     l = (53 * l + 1) % 169;
  460.     if ((l * m) % 64 >= 32)
  461.        s += t;
  462.     t *= 0.5;
  463.       }
  464.       x[ii] = (long) floor(s * NO_NUMBERS);
  465.    }
  466.  
  467.    c = 362436;
  468.  
  469.    k_ptr = &x[FIRSTDIM - 1];
  470.    j_ptr = &x[LASTDIM - 1];
  471. }
  472. #endif
  473.  
  474.